Optimization of a State-to-State Transfer in a Two-Level-System¶
In [1]:
# NBVAL_IGNORE_OUTPUT
%load_ext watermark
#%load_ext autoreload
#%autoreload 2
import qutip
import numpy as np
import scipy
import matplotlib
import matplotlib.pylab as plt
import krotov
import os
from scipy.signal import savgol_filter
%matplotlib notebook
%watermark -v --iversions
π = np.pi
sqrt = np.sqrt
basis = qutip.basis
tensor = qutip.tensor
coherent = qutip.coherent
$\newcommand{tr}[0]{\operatorname{tr}} \newcommand{diag}[0]{\operatorname{diag}} \newcommand{abs}[0]{\operatorname{abs}} \newcommand{pop}[0]{\operatorname{pop}} \newcommand{aux}[0]{\text{aux}} \newcommand{opt}[0]{\text{opt}} \newcommand{tgt}[0]{\text{tgt}} \newcommand{init}[0]{\text{init}} \newcommand{lab}[0]{\text{lab}} \newcommand{rwa}[0]{\text{rwa}} \newcommand{bra}[1]{\langle#1\vert} \newcommand{ket}[1]{\vert#1\rangle} \newcommand{Bra}[1]{\left\langle#1\right\vert} \newcommand{Ket}[1]{\left\vert#1\right\rangle} \newcommand{Braket}[2]{\left\langle #1\vphantom{#2} \mid
2\vphantom{#1}\right\rangle}¶
\newcommand{op}[1]{\hat{#1}} \newcommand{Op}[1]{\hat{#1}} \newcommand{dd}[0]{\,\text{d}} \newcommand{Liouville}[0]{\mathcal{L}} \newcommand{DynMap}[0]{\mathcal{E}} \newcommand{identity}[0]{\mathbf{1}} \newcommand{Norm}[1]{\lVert#1\rVert} \newcommand{Abs}[1]{\left\vert#1\right\vert} \newcommand{avg}[1]{\langle#1\rangle} \newcommand{Avg}[1]{\left\langle#1\right\rangle} \newcommand{AbsSq}[1]{\left\vert#1\right\vert^2} \newcommand{Re}[0]{\operatorname{Re}} \newcommand{Im}[0]{\operatorname{Im}}$ The purpose of this example is not to solve an especially interesting physical problem but to give a rather simple example of how the package can be used in order to solve an optimization problem.
Define the Hamiltonian¶
In the following the Hamiltonian, guess field and states are defined.
The Hamiltonian $\op{H}_{0} = - \omega \op{\sigma}_{z}$ represents a simple qubit with energy level splitting $\omega$ in the basis $\{\ket{0},\ket{1}\}$. The control field $\epsilon(t)$ is assumed to couple via the Hamiltonian $\op{H}_{1}(t) = \epsilon(t) \op{\sigma}_{x}$ to the qubit, i.e., the control field effectively drives transitions between both qubit states. For now, we initialize the control field as constant.
In [397]:
L = 3
N = 10
α = 2
Plotting functions¶
In [475]:
def to_two_level(state):
if state.type is 'oper':
return qutip.Qobj(state[0:2,0:2])
else:
return qutip.Qobj(state[0:2])
def plot_population(n, tlist):
fig, ax = plt.subplots(figsize=(15,4))
leg = []
for i in range(len(n)):
ax.plot(tlist, n[i], label=str(i))
leg.append(str(i))
ax.legend()
ax.set_xlabel('Time (ns)')
ax.set_ylabel('Occupation')
ax.legend(leg)
plt.show(fig)
def plot_pulse(pulse, tlist, T=None, ylim=False):
fig, ax = plt.subplots(figsize=(15,4))
if callable(pulse):
pulse = np.array([pulse(t, args=None) for t in tlist])
if np.any(np.iscomplex(pulse)):
ax.plot(tlist, np.real(pulse))
ax.plot(tlist, np.imag(pulse))
ax.legend(['Re', 'Im'])
else:
ax.plot(tlist, pulse)
if T is not None:
ax.plot(tlist, [S(t, T) for t in tlist], color='k', linestyle='--', linewidth=1)
ax.plot(tlist, [-S(t, T) for t in tlist], color='k', linestyle='--', linewidth=1)
ax.set_xlabel('Time (ns)')
ax.set_ylabel('Pulse amplitude')
if ylim is False:
ax.set_ylim([-pulse_max(0)*1.05,pulse_max(0)*1.05])
plt.show(fig)
def plot_system(ψ):
bl = qutip.Bloch()
bl.add_states(to_two_level(ψ.ptrace(0)))
bl.show()
qutip.visualization.plot_wigner_fock_distribution(to_two_level(ψ.ptrace(1)))
def plot_resonator(ψ):
fig, ax = plt.subplots(1,len(ψ), figsize=(3*len(ψ),3))
if len(ψ)<=1:
ψ = ψ[0]
axis = ax
qutip.visualization.plot_wigner(ψ.ptrace(1), fig=fig, ax=axis, alpha_max = 2*α)
axis.axis_equal = True
else:
for (ϕ, axis) in zip(ψ, ax):
qutip.visualization.plot_wigner(ϕ.ptrace(1), fig=fig, ax=axis, alpha_max = 2*α)
axis.axis_equal = True
def plot_cardinal(ψ):
bl = qutip.Bloch()
bl.vector_color = ['r','g','b','g','b','r']
[bl.add_states(to_two_level(ϕ.ptrace(0)), 'vector') for ϕ in ψ]
bl.show()
def plot_all(dyn, ψ):
ψ_i = [g.states[0] for g in dyn]
ψ_f = [g.states[-1] for g in dyn]
ψ_t = [ϕ[1] for ϕ in ψ]
plot_cardinal(ψ_i)
plot_resonator(ψ_i)
plot_cardinal(ψ_t)
plot_resonator(ψ_t)
plot_cardinal(ψ_f)
plot_resonator(ψ_f)
def plot_evolution(dyn, steps=1):
for d in dyn:
points = [to_two_level(s.ptrace(0)) for s in d.states[0:-1:steps]]
bl = qutip.Bloch()
bl.vector_color = 'r'
bl.point_color = 'r'
bl.point_marker = 'o'
bl.add_states(points, 'point')
bl.show()
bl = qutip.Bloch()
bl.vector_color = 'r'
bl.point_color = 'r'
bl.point_marker = 'o'
bl.view = [bl.view[0], 80]
bl.add_states(points, 'point')
bl.show()
def get_objectives(T=None):
#if use_rotating:
# objectives = [krotov.Objective(initial_state=ψ[0], target=ψ[1], H=H) for ψ in state_rot(ϕ, T)]
#else:
objectives = [krotov.Objective(initial_state=ψ[0], target=ψ[1], H=H) for ψ in ϕ]
return objectives
def plot_matrix_final_target(target_state, final_state, xlabels, ylabels, el=30, az=135):
fig, ax = qutip.visualization.matrix_histogram(final_state * target_state.dag(), xlabels, ylabels, colorbar=False, limits=[-1,1])
qutip.visualization.matrix_histogram(proj(target_state), xlabels, ylabels, colorbar=False, limits=[-1,1], fig=fig, ax=ax)
facecolors = np.zeros((6*L**2,4))*0.1
edgecolors = np.tile([0,0,0,0.9], (6*L**2,1))
ax.get_children()[2].set_facecolors(facecolors)
ax.get_children()[2].set_edgecolors(edgecolors)
ax.set_zticks(np.arange(-1,1,0.25))
ax.view_init(elev=el, azim=az)
return (fig, ax)
In [399]:
def fid(result, target):
return (np.abs((result.states[-1].dag()*target).full())**2)[0][0]
def proj(ψ, ϕ=None):
if ϕ is None:
return ψ * ψ.dag()
else:
return ψ * ϕ.dag()
In [411]:
cat_0 = (coherent(N, α) + coherent(N,-α)).unit()
cat_1 = (coherent(N, 1j*α) + coherent(N,-1j*α)).unit()
Si = qutip.operators.identity(L)
Ri = qutip.operators.identity(N)
I = tensor(Si,Ri)
zero_q = qutip.operators.qzero(L)
zero_r = qutip.operators.qzero(N)
zero_q = tensor(zero_q, Ri)
zero_r = tensor(Si, zero_r)
zero = zero_q + zero_r
σ_z = proj(qutip.basis(L, 0)) - proj(qutip.basis(L, 1))
σ_y = 1j*(proj(qutip.basis(L, 1),qutip.basis(L, 0)) - proj(qutip.basis(L, 0), qutip.basis(L, 1)))
σ_x = proj(qutip.basis(L, 0),qutip.basis(L, 1)) - proj(qutip.basis(L, 1), qutip.basis(L, 0))
σ_z = tensor(σ_z, Ri)
σ_y = tensor(σ_y, Ri)
σ_x = tensor(σ_x, Ri)
b = qutip.operators.destroy(L)
a = qutip.operators.destroy(N)
b = tensor(b, Ri)
a = tensor(Si, a)
ω_r = 8.3056 * 2 * π # resonator frequency
ω_q = 6.2815 * 2 * π # qubit frequency
K_q = -2*π*297e-3 # Kerr qubit 200-300 MHz
K_r = 2*π*4.5e-6 # Kerr res 1-10 Khz
ω_ef = ω_q + K_q
ω_gf = ω_q + K_q/2
χ = 0.025 * 2 * π # parameter in the dispersive hamiltonian
Δ = abs(ω_r - ω_q) # detuning
g = sqrt(Δ * χ) # coupling strength that is consistent with chi
γ = 1e1
use_rotating = True
def hamiltonian(ω=1.0, ampl0=1, use_rotating=True, pulses=None, tlist=None, start_pulse=None, T=1, phase=0, σ=σ):
"""Two-level-system Hamiltonian
Args:
ω (float): energy separation of the qubit levels
ampl0 (float): constant amplitude of the driving field
"""
#H_occ = w_r*a.dag()*a + w_q*b.dag()*b
if L==3:
H_occ_q = qutip.Qobj(np.diag([0, ω_q, 2*ω_q]))
else:
H_occ_q = qutip.Qobj(np.diag([0, ω_q]))
H_occ_q = tensor(H_occ_q, Ri)
H_occ_r = ω_r * a.dag()*a
H_occ = H_occ_q + H_occ_r
use_dispersive = True
use_kerr = True
#if use_dispersive:
# H_coup = - χ * a.dag()*a * b.dag()*b
#H_coup = χ * (a.dag()*a + I/2) * σ_z
#else:
# H_coup = g * (a.dag() * b + a * b.dag())
# H_coup = g * σ_x *a.dag() + a
if use_kerr:
H_kq = - K_q/2 * b.dag()**2 * b**2
H_kr = - K_r/2 * a.dag()**2 * a**2
else:
H_kq = zero_q
H_kr = zero_r
H_d = zero_q + zero_q
if use_rotating:
H_d += H_kq + H_kr #+ H_coup
H_qr = (b.dag() + b)
H_qi = 1j*(b.dag() - b)
H_rr = (a.dag() + a)
H_ri = 1j*(a.dag() - a)
H_cc = g*(a*b + a.dag()*b.dag())
H_cs = g*1j*(a*b - a.dag()*b.dag())
if start_pulse is None:
ϵ_qr = lambda t, args: ampl0
ϵ_qi = lambda t, args: ampl0
ϵ_rr = lambda t, args: ampl0
ϵ_ri = lambda t, args: ampl0
else:
ϵ_qr = shape_field(lambda t, args: ampl0, start_pulse, T, σ)
ϵ_qi = shape_field(lambda t, args: ampl0, start_pulse, T, σ)
ϵ_rr = shape_field(lambda t, args: ampl0, start_pulse, T, σ)
ϵ_ri = shape_field(lambda t, args: ampl0, start_pulse, T, σ)
if pulses:
ϵ_qr = pulses[0]
ϵ_qi = pulses[1]
ϵ_rr = pulses[2]
ϵ_ri = pulses[3]
ϵ_cc = lambda t, args: np.cos((ω_q - ω_r)*t)
ϵ_cs = lambda t, args: np.sin((ω_q - ω_r)*t)
return [H_d, [H_qr, ϵ_qr], [H_qi, ϵ_qi], [H_rr, ϵ_rr], [H_ri, ϵ_ri], [H_cc, ϵ_cc], [H_cs, ϵ_cs]]
else:
H_d += H_occ + H_kerr#+ H_coup
H_q = b
H_qc = b.dag()
H_rr = ZERO
H_ri = ZERO
ϵ_q = lambda t, args: 1j*ampl0*np.exp(1j*ω_q*t)
ϵ_qc = lambda t, args: -1j*ampl0*np.exp(-1j*ω_q*t)
#ϵ_rr = lambda t, args: ampl0
#ϵ_ri = lambda t, args: ampl0
if pulses:
ϵ_q = pulses[0]
ϵ_qc = pulses[1]
#ϵ_rr = np.zeros(len(pulses[0]))
#ϵ_ri = np.zeros(len(pulses[0]))
return [H_d, [H_q, ϵ_q], [H_qc, ϵ_qc]]#, [H_rr, ϵ_rr], [H_ri, ϵ_ri]]
def coeffs_to_state(c,init = True):
if init:
ψ = tensor((c[0]*basis(L,0) + c[1]*basis(L,1)).unit() , (basis(N,0)))
else:
ψ = tensor((basis(L,0)) , (c[0]*cat_0 + c[1]*cat_1).unit())
return ψ
def states(coeffs):
return [[coeffs_to_state(c,True),coeffs_to_state(c,False)] for c in coeffs]
Define the optimization target¶
First we define the time grid of the dynamics, i.e., by taking the following values as an example, we define the initial state to be at time $t=0$ and consider a total propagation time of $T=5$. The entire time grid is divided into $n_{t}=500$ equidistant time steps.
In [412]:
def pulse_max(σ):
A = 1.56246130414 # Chosen such that the integral of any Blackman pulse = π
A = A/2
σ = np.max((σ,3))
return A/(np.sqrt(2*π)*σ)
In [413]:
σ_max = 3 # ns (gaussian pulse limit)
amp_max = pulse_max(0)
T_q = (2*π)/ω_q
T = 18
σ = T/6
steps = 4**3*int(np.ceil(T))
tlist = np.linspace(0, T, steps)
H = hamiltonian(ampl0=1, use_rotating=True, phase=np.exp(-1j*ω_q*T))
#coeffs = [(1,0), (1,-1), (1,1j), (1,1), (1,-1j), (0,1)]
coeffs = [(1,0)]
st = states(coeffs)
ϕ = [[ tensor(basis(L,1), basis(N,0)), tensor((basis(L,0)).unit(), (coherent(N,1))) ]]
F_err = 1e-5
F_oc_tar = 1-F_err
Next, we define the optimization targets, which is technically a list of objectives, but here it has just one entry defining a simple state-to-state transfer from initial state $\ket{\Psi_{\init}} = \ket{0}$ to the target state $\ket{\Psi_{\tgt}} = \ket{1}$, which we want to reach at final time $T$. Note that we also have to pass the Hamiltonian $\op{H}(t)$ that determines the dynamics of the system to the optimization objective.
In [414]:
import copy
def state_rot(ϕ, T):
ϕ = copy.deepcopy(ϕ)
ϕ_q = ϕ.ptrace(0)
if np.sum(np.array(ϕ[0][1].full())==0) != L-1:
if L == 3:
rot_evo = qutip.Qobj([[1, 0, 0],[0, np.exp(-1j * ω_q * T), 0],[0, 0, 0]])
else:
rot_evo = qutip.Qobj([[1, 0],[0, np.exp(-1j * ω_q * T)]])
ϕ[0][1] = rot_evo * ϕ[0][1]
#ϕ_r = ϕ.ptrace(0)
return ϕ
if use_rotating:
objectives = [krotov.Objective(initial_state=ψ[0], target=ψ[1], H=H) for ψ in ϕ]
#objectives = [krotov.Objective(initial_state=ψ[0], target=ψ[1], H=H) for ψ in state_rot(copy.deepcopy(ϕ), T)]
else:
objectives = [krotov.Objective(initial_state=ψ[0], target=ψ[1], H=H) for ψ in ϕ]
In addition, we have to define and assign a shape function $S(t)$ for the update in each control iteration to each control field that will be updated. This shape usually takes care of experimental limits such as the necessity of finite ramps at the beginning and end of the control field or other conceivable limitations for field shapes: wherever $S(t)$ is zero, the optimization will not change the value of the control from the original guess.
At this point, we also change the initial control field $\epsilon_{0}(t)$ from a constant to a shaped pulse that switches on smoothly from zero and again switches off at the final time $T$. We re-use the shape function $S(t)$ that we defined for the updates for this purpose (although generally, $S(t)$ for the updates has nothing to with the shape of the control field).
In [415]:
def S(t, T=6*σ, σ=σ):
rise_time = 2
return amp_max*krotov.shapes.flattop(t, t_start=0, t_stop=T, t_rise=rise_time, t_fall=rise_time, func='sinsq')
def start_pulse(t, T=6*σ, σ=σ):
if σ is not None:
return pulse_max(σ)*krotov.shapes.blackman(t, t_start=0, t_stop=6*σ)
else:
return krotov.shapes.blackman(t, t_start=0, t_stop=T)
def zero_pulse(t, T=T, σ=4):
return 0
def unit_pulse(t, T=T, σ=4):
return 1
def shape_field(ϵ, sf, T, σ):
"""Applies the shape function S(t) to the guess field"""
ϵ_shaped = lambda t, args: ϵ(t, args)*sf(t, T=T, σ=σ)
return ϵ_shaped
S_unit = [unit_pulse]*4
S_zero = [zero_pulse]*4
#S_start = [lambda t, T=T, σ=σ: 0.01*unit_pulse(t, T=T, σ=σ),start_pulse]
S_start = [zero_pulse, start_pulse, zero_pulse, start_pulse]
S_funs = [S]*4
for i, H_i in enumerate(H[1:5]):
H_i[1] = shape_field(H_i[1], S_start[i], T, σ)
Having defined the shape function $S(t)$ and having shaped the guess field, we
now tell the optimization to also use $S(t)$ as the update-shape for
$\epsilon_0(t)$. In addition, we have to choose lambda_a for each control
field. It controls the update magnitude of the respective field in each
iteration.
It is convenient to introduce the function print_fidelity, which can be passed
to the optimization procedure and will be called after each iteration and thus
provides additional feedback about the optimization progress.
Simulate dynamics of the guess field¶
Before heading towards the optimization procedure, we first simulate the dynamics under the guess field $\epsilon_{0}(t)$.
The following plot shows the guess field $\epsilon_{0}(t)$, which is, as chosen above, just a constant field (with a smooth switch-on and switch-off)
In [416]:
for i, H_i in enumerate(H[1:]):
plot_pulse(H_i[1], tlist, ylim=(i > 3))
The next line solves the equation of motion for the defined objective, which contains the initial state $\ket{\Psi_{\init}}$ and the Hamiltonian $\op{H}(t)$ defining its evolution.
In [417]:
guess_dynamics = [ob.mesolve(tlist, progress_bar=True, options=qutip.Options(nsteps=50000)) for ob in objectives]
In [418]:
def qubit_occupation(dyn):
occ = [basis(L,i)*basis(L,i).dag() for i in range(0,L)]
states = [state.ptrace(0) for state in dyn.states]
n = qutip.expect(occ, states)
plot_population(n, dyn.times)
def resonator_occupation(dyn):
occ = [basis(N,i)*basis(N,i).dag() for i in range(0,N)]
states = [state.ptrace(1) for state in dyn.states]
n = qutip.expect(occ, states)
plot_population(n, dyn.times)
def plot_norm(result):
state_norm = lambda i: result.states[i].norm()
states_norm=np.vectorize(state_norm)
fig, ax = plt.subplots()
ax.plot(result.times, states_norm(np.arange(len(result.states))))
ax.set_title('Norm loss', fontsize = 15)
ax.set_xlabel('Time (ns)')
ax.set_ylabel('State norm')
plt.show(fig)
In [419]:
def desuperposition(ϕ, f):
ϕ_f = copy.deepcopy(ϕ)
for (i, s) in enumerate(ϕ_f.full()):
if np.abs(s)[0] < f:
ϕ_f.data[i] = 0
ϕ_f = ϕ_f.unit()
for (i, s) in enumerate(ϕ_f.full()):
if np.abs(s)[0] == 1:
ϕ_f.data[i] = 1
return ϕ_f
In [420]:
final_state = guess_dynamics[0].states[-1]
dm = final_state * ϕ[0][1].dag()
fid = np.abs((final_state.dag() * ϕ[0][1]).full()[0][0])**2
print('Fidelity: {}'.format(fid))
#qutip.visualization.matrix_histogram_complex(dm, limits=[0,1])
plot_cardinal([ϕ[0][0].ptrace(0), ϕ[0][1].ptrace(0), final_state.ptrace(0)])
qutip.visualization.plot_wigner_fock_distribution(ϕ[0][1].ptrace(1))
qutip.visualization.plot_wigner_fock_distribution(final_state.ptrace(1))
Out[420]:
In [421]:
qubit_occupation(guess_dynamics[0])
resonator_occupation(guess_dynamics[0])
#plot_norm(guess_dynamics[0])
The plot of the population dynamics shows that the guess field does not transfer the initial state $\ket{\Psi_{\init}} = \ket{0}$ to the desired target state $\ket{\Psi_{\tgt}} = \ket{1}$.
In [422]:
#plot_all(guess_dynamics, ϕ)
plot_evolution(guess_dynamics, steps=5)
In [423]:
qubit_pulses = [H[2][1](t, 0) for t in tlist]
res_pulses = [H[4][1](t, 0) for t in tlist]
#qubit_pulses_filtered = apply_spectral_filter(copy.deepcopy(qubit_pulses), tlist, 0, 0.5)
plot_spectrum(qubit_pulses, tlist, mark_freq=[0, -K_q, -K_q/2], pos=0, xlim=[-2*π,2*π], title="Qubit")
plot_spectrum(res_pulses, tlist, mark_freq=[0, -K_r, -K_r/2], pos=0, xlim=[-π/4,π/4], title="Resonator")
#plot_spectrum(qubit_pulses_filtered, tlist, mark_freq=[0, K_q], pos=0, xlim=[-2*K_q,2*K_q])
#plot_pulse(qubit_pulses, tlist)
#plot_pulse(qubit_pulses_filtered, tlist)
#plot_spectrum(qubit_pulses[1], time_list, mark_freq=[ω_q, ω_q + K_q, ω_q - K_q], pos=0)
Optimize¶
In the following we optimize the guess field $\epsilon_{0}(t)$ such that the intended state-to-state transfer $\ket{\Psi_{\init}} \rightarrow \ket{\Psi_{\tgt}}$ is solved.
The cell below carries out the optimization. It
requires, besides the
previously
defined optimization objectives, information
about the
optimization functional
$F$ (via chi_constructor) and the
propagation method that should be used. In
addition, the number of total
iterations is required and, as an option, we pass
an info-hook that after each
iteration combines a complete printout of the state
of the optimization with the
print_fidelity function defined above.
Here, we choose $F = F_{re}$ with \begin{equation}
F_{re}¶
\Re\Braket{\Psi(T)}{\Psi_{\tgt}} \end{equation}
with $\ket{\Psi(T)}$ the forward propagated state of $\ket{\Psi_{\init}}$.
In [97]:
pulse_options = {H_i[1]: dict(lambda_a=5, shape=S_funs[i]) for i, H_i in enumerate(H[1:])}
In [98]:
ω_0 = 0
ω_1 = np.abs(K_q/2)
In [478]:
def F_oc(fw_states_T, objectives, tau_vals=None, **kwargs):
return np.abs(krotov.functionals.f_tau(fw_states_T, objectives, tau_vals, **kwargs))**2
def calc_fidelity(tau_vals):
return np.abs(np.sum(tau_vals)/len(tau_vals))**2
def print_fidelity(**args):
fid = calc_fidelity(np.array(args['tau_vals']))
print(" F_t = {} | F = {} | F_t - F = {}".format(F_oc_tar, fid, F_oc_tar-fid))
def plot_fid_convergence(info_vals):
fig, ax = plt.subplots(1,1)
ax.plot(info_vals)
ax.set_xticks(np.arange(0, len(info_vals), step=1))
ax.set_xlabel('Iteration')
ax.set_ylabel('Fidelity')
#ax.set_ylim((-0.2,.2))
plt.show()
def plot_fid_convergence(ax, info_vals, T):
ax.plot3D(range(0,len(info_vals)), [T]*len(info_vals), info_vals)
def modify_params(**kwargs):
spectral_filtering = False
if spectral_filtering:
# Spectral filtering
tlist = kwargs['tlist']
dt = tlist[1] - tlist[0] # assume equi-distant time grid
n = len(tlist) - 1 # = len(pulse)
# remember that pulses are defined on intervals of tlist
ω = np.abs(np.fft.fftfreq(n, d=dt / (2.0 * π)))
# the normalization factor 2π means that w0 and w1 are angular
# frequencies, corresponding directly to energies in the Hamiltonian
# (ħ = 1).
flt = (ω_0 <= ω) * (ω <= ω_1)
# flt is the (boolean) filter array, equivalent to an array of values 0
# and 1
shape_arrays = kwargs['shape_arrays']
for (i, (pulse, shape)) in enumerate(zip(kwargs['optimized_pulses'], shape_arrays)):
spectrum = np.fft.fft(pulse)
# apply the filter by element-wise multiplication
spectrum[:] *= flt[:]
# after the inverse fft, we should also multiply with the
# update shape function. Otherwise, there is no guarantee that
# the filtered pulse will be zero at t=0 and t=T (assuming that
# is what the update shape is supposed to enforce). Also, it is
# important that we overwrite `pulse` in-place (pulse[:] = ...)
kwargs['optimized_pulses'][i][:] = np.fft.ifft(spectrum).real * shape
# Limit pulse amplitude to 1 and smoothen
for i in range(len(kwargs['optimized_pulses'])-2):
#kwargs['optimized_pulses'][i] = savgol_filter(kwargs['optimized_pulses'][i], 9, 2)
pulse_max = np.max(np.abs(kwargs['optimized_pulses'][i]))
if pulse_max > amp_max:
kwargs['optimized_pulses'][i] = (amp_max*np.array(kwargs['optimized_pulses'][i])/pulse_max)
kwargs['optimized_pulses'][i] = np.fmax(np.fmin(kwargs['optimized_pulses'][i], kwargs['shape_arrays'][i]), -np.array(kwargs['shape_arrays'][i]))
#conv = 3*σ
#if (conv % 2 == 0): conv += 1
#kwargs['optimized_pulses'][i] = savgol_filter(kwargs['optimized_pulses'][i], conv, 2)
for i in range(len(kwargs['optimized_pulses'])):
if kwargs['iteration'] % 10 == 0:
plot_pulse(kwargs['optimized_pulses'][i], kwargs['tlist'][:-1], kwargs['tlist'][-1])
#plot_spectrum(kwargs['optimized_pulses'][i], kwargs['tlist'][:-1], mark_freq=[0, -K_q, -K_q/2], mark_color=['r','g','b'], pos=0, xlim=[-(2*π), (2*π)])
# Update λ
#fac = 1
#steps = 5
#λₐ = kwargs['lambda_vals'][0]
#for i in range(len(kwargs['lambda_vals'])):
# kwargs['lambda_vals'][i] = λₐ * fac
# lambda_a = λₐ * fac
#print("λₐ = {}".format(kwargs['lambda_vals']))
In [479]:
# Reset results
opt_result = None
In [480]:
def convergence_reason(opt_result):
if opt_result == None:
return True
reasons = ['monotonic', 'iterations']
for r in reasons:
if opt_result.message.find(r)>0:
return r
if opt_result.message.find('F_oc')>0 or opt_result.message.find('Δ')>0:
return False
In [481]:
from datetime import datetime
def current_time():
return datetime.now().isoformat()[:16].replace('T',' ')
In [482]:
def sine_pulse(t, T=σ*6, σ=σ):
ω = K_q/2
return np.sin(ω*t)*start_pulse(t, T=T, σ=σ)
In [483]:
def run_optim(T, lambda_a, ϕ):
σ = T/6
opt_result = None
print('T = {}'.format(T))
#steps = 500
tlist = np.linspace(0, T, 16*int(np.ceil(T)))
s_pulse = None
H = hamiltonian(ampl0=1, use_rotating=True, start_pulse=s_pulse, T=T, phase=np.exp(-1j*ω_q*T))
S_start = [start_pulse]*4
S_funs = [S]*4
for i, H_i in enumerate(H[1:]):
if i < 4:
H_i[1] = shape_field(H_i[1], S_start[i], T, σ)
plot_pulse(H_i[1], tlist)
else:
plot_pulse(H_i[1], tlist, ylim=True)
#objectives = [krotov.Objective(initial_state=ψ[0], target=ψ[1], H=H) for ψ in state_rot(ϕ, T)]
objectives = [krotov.Objective(initial_state=ψ[0], target=ψ[1], H=H) for ψ in ϕ]
# Check if guess pulse realises
guess_dynamics = [ob.mesolve(tlist, options=qutip.Options(nsteps=50000)) for ob in objectives]
final_state = guess_dynamics[0].states[-1]
dm = final_state * ϕ[0][1].dag()
fid = np.abs((final_state.dag() * ϕ[0][1]).full()[0][0])**2
if fid > F_oc_tar:
print('Guess pulse realises transfer already.')
return True
pulse_options = {H_i[1]: dict(lambda_a=lambda_a, shape=lambda t: S_funs[i](t, T=T, σ=σ)) for i, H_i in enumerate(H[1:5])}
pulse_options.update({H_i[1]: dict(lambda_a=lambda_a, shape=lambda t: 0) for i, H_i in enumerate(H[5:])})
#pulse_options = {
# H[2][1]: dict(lambda_a=lambda_a, shape=0),
# H[1][1]: dict(lambda_a=lambda_a, shape=lambda t: S_funs[0](t, T=T, σ=σ)),
#}
while convergence_reason(opt_result):
if convergence_reason(opt_result) == 'monotonic':
break
#lambda_a *= 2
# print('λₐ = {}'.format(lambda_a))
# pulse_options = {H_i[1]: dict(lambda_a=lambda_a, shape=lambda t: S_funs[i](t, T)) for i, H_i in enumerate(H[1:])}
#iters = 5
#if opt_result is not None:
# iters = opt_result.iters[-1] + iters
opt_result = krotov.optimize_pulses(
objectives,
pulse_options=pulse_options,
tlist=tlist,
propagator=krotov.propagators.expm,
chi_constructor=krotov.functionals.chis_ss,
info_hook=krotov.info_hooks.chain(
krotov.info_hooks.print_table(J_T=F_oc),
print_fidelity
),
check_convergence=krotov.convergence.Or(
krotov.convergence.value_above(F_oc_tar, name='F_oc'),
krotov.convergence.delta_below(F_err*1e-1),
#krotov.convergence.check_monotonic_fidelity,
),
modify_params_after_iter = modify_params,
#iter_stop=1,
continue_from = opt_result,
)
print(opt_result.message)
opt_result.dump(os.path.join(os.getcwd(),'results','{}_pi_pulse_optim_{}.dat'.format(current_time(),T)))
In [484]:
step_size = pulse_max(0)*2.5
λ = 1/step_size
ϕ = [[ tensor(basis(L,1), basis(N,0)), tensor((basis(L,0)).unit(), cat_1) ]]
existing_times = [float(file.split('_')[4][:-4]) for file in os.listdir('results')]
t_times = np.flip(np.arange(1,21.5,1))
t_times = [42.]
for tot in t_times:
if tot not in [float(file.split('_')[4][:-4]) for file in os.listdir('results')]:
#plot_cardinal(state_rot(ϕ, tot)[0])
if tot.is_integer():
tot = int(tot)
run_optim(tot, λ, ϕ)
else:
print('T = {} already exists'.format(tot))
Simulate dynamics of the optimized field¶
Having obtained the optimized control field, we can now plot it and calculate the population dynamics under this field.
In [465]:
def plot_results_3d(results):
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.set_xlabel('Iteration')
ax.set_zlabel('Fidelity')
ax.set_ylabel('Pulse length')
ax.set_zlim(0,1.1)
for (r, T) in results:
plot_fid_convergence(ax, r.info_vals[1:], T)
ax.view_init(elev=20, azim=340)
return (fig, ax)
def plot_results_iteration(results):
fig = plt.figure()
ax = plt.axes()
ax.set_xlabel('Iteration')
ax.set_ylabel('Fidelity')
for (r, T) in results:
ax.plot(range(0,len(r.info_vals)-1), r.info_vals[1:])
#print('F = {}'.format(r.info_vals[-1]))
return (fig, ax)
def plot_results_pulse_length(results, iteration=-1, ax=None):
if ax is None:
ax = plt.axes()
else:
ax.clear()
ax.set_xlabel('Pulse length')
ax.set_ylabel('Fidelity')
for (r, T) in results:
it = min(len(r.info_vals)-1,iteration)
ax.plot(T, r.info_vals[it], 'ok')
ax.set_ylim(0,1.1)
return ax
In [466]:
folder = 'results'
results = [(krotov.result.Result.load(os.path.join(os.getcwd(),folder,file), objectives=get_objectives(T=float(file.split('_')[-1][:-4]))), float(file.split('_')[-1][:-4])) for file in os.listdir(folder) if file[-4:]=='.dat']
%matplotlib inline
plot_results_3d(results)
#plot_results_pulse_length(results, iteration=0)
plot_results_iteration(results)
Out[466]:
In [565]:
plot_results_pulse_length(results, iteration=20, ax=ax)
ax = plt.axes()
def interactive_plot(iteration):
plot_results_pulse_length(results, iteration=iteration, ax=ax)
interact(interactive_plot, iteration=widgets.IntSlider(min=0,max=900,step=1,value=0));
In [1394]:
import matplotlib2tikz
matplotlib2tikz.save("mytikz.tex")
In [1486]:
%matplotlib inline
Analyze¶
In [206]:
results = [(krotov.result.Result.load(os.path.join(os.getcwd(),folder,file), objectives=get_objectives(T=float(file.split('_')[-1][:-4]))), float(file.split('_')[-1][:-4])) for file in os.listdir(folder) if file[-4:]=='.dat']
results = results[-1:]
In [207]:
steps2 = len(results[0][0].tlist)*1000
for (r,_) in results:
tlist = r.tlist
#opt_dynamics = [ob.mesolve(tlist, progress_bar=True) for ob in r.objectives]
#qubit_occupation(opt_dynamics[0])
c = r.optimized_controls
tlist2 = np.linspace(0, tlist[-1], steps2)
Ω_q = c[0]+1j*c[1]
Ω_q = np.interp(tlist2, tlist, Ω_q)
Ω_r = c[2]+1j*c[3]
Ω_r = np.interp(tlist2, tlist, Ω_r)
pulses_lab = [Ω_q*np.exp(1j*ω_q*tlist2), np.conj(Ω_q)*np.exp(-1j*ω_q*tlist2), Ω_r*np.exp(1j*ω_r*tlist2), np.conj(Ω_r)*np.exp(-1j*ω_r*tlist2)]
opt_dynamics = [ob.mesolve(tlist, progress_bar=True) for ob in r.optimized_objectives]
plot_pulse(r.guess_controls[0], tlist)
plot_pulse(r.guess_controls[1], tlist)
plot_pulse(r.guess_controls[2], tlist)
plot_pulse(r.guess_controls[3], tlist)
plot_pulse(c[0], tlist)
plot_pulse(c[1], tlist)
plot_pulse(c[2], tlist)
plot_pulse(c[3], tlist)
plot_pulse(pulses_lab[0], tlist2)
plot_pulse(pulses_lab[1], tlist2)
plot_pulse(pulses_lab[2], tlist2)
plot_pulse(pulses_lab[3], tlist2)
qubit_occupation(opt_dynamics[0])
resonator_occupation(opt_dynamics[0])
plot_spectrum(pulses_lab[0], tlist2, mark_freq=[ω_q, ω_ef, ω_gf],mark_color=['r','g','b'], pos=0, xlim=[ω_q*0.9, ω_q*1.1])
plot_spectrum(pulses_lab[2], tlist2, mark_freq=[ω_r],mark_color=['r','g','b'], pos=0, xlim=[ω_r*0.9, ω_r*1.1])
#plot_spectrum(pulses_lab[1], tlist2, mark_freq=[ω_q, ω_ef, ω_gf], pos=0, xlim=[-ω_q*0.95, -ω_q*1.05])
#H_lab = hamiltonian(ampl0=1, use_rotating=False, pulses=pulses_lab)
#objectives_lab = [krotov.Objective(initial_state=ψ[0], target=ψ[1], H=H_lab) for ψ in ϕ]
In [720]:
xlabels = ['$|0\\rangle$','$|1\\rangle$','$|2\\rangle$']
ylabels = ['$\\langle 0|$','$\\langle 1|$','$\\langle 2|$']
final_state = desuperposition(opt_dynamics[0].states[-1], F_err)
#final_state = opt_dynamics[0].states[-1]
target_state = results[0][0].objectives[0].target
plot_matrix_final_target(target_state, final_state, xlabels, ylabels, el=45, az=150)
plot_matrix_final_target(target_state, final_state, xlabels, ylabels, el=10, az=150)
plot_cardinal([target_state, final_state])
plot_evolution(opt_dynamics)
Something something¶
In [645]:
[plot_pulse(c, tlist) for c in opt_result.optimized_controls];
In [ ]:
plot_pulse(np.abs(opt_result.optimized_controls[0]+1j*opt_result.optimized_controls[1]), tlist)
plot_pulse(np.angle(opt_result.optimized_controls[0]+1j*opt_result.optimized_controls[1]), tlist)
In contrast to the dynamics under the guess field, the optimized field indeed drives the initial state $\ket{\Psi_{\init}} = \ket{0}$ to the desired target state $\ket{\Psi_{\tgt}} = \ket{1}$.
In [ ]:
opt_dynamics = [ob.mesolve(tlist, progress_bar=True) for ob in opt_result.optimized_objectives]
In [ ]:
qubit_occupation(opt_dynamics[0])
In [ ]:
plot_evolution(opt_dynamics)
#plot_all(opt_dynamics, ϕ)
Simulate dynamics in lab frame¶
In [ ]:
#rotating_pulses = [np.vectorize(H_i[1])(tlist, None) for H_i in H[1:]]
rotating_pulses = oct_result.optimized_controls
tlist2 = np.linspace(0, T, int(np.ceil(200*T/T_q)))
Ω = rotating_pulses[0]+1j*rotating_pulses[1]
Ω = np.interp(tlist2, tlist, Ω)
pulses_lab = [Ω*np.exp(1j*ω_q*tlist2), np.conj(Ω)*np.exp(-1j*ω_q*tlist2)]
H_lab = hamiltonian(ω=1.0, ampl0=1, use_rotating=False, pulses=pulses_lab)
objectives_lab = [krotov.Objective(initial_state=ψ[0], target=ψ[1], H=H_lab) for ψ in ϕ]
opt_dynamics_lab = [ob.mesolve(tlist2, progress_bar=True, options=qutip.Options(nsteps=50000)) for ob in objectives_lab]
In [ ]:
plot_pulse(Ω, tlist2)
[plot_pulse(p, tlist2) for p in pulses_lab];
In [ ]:
qubit_occupation(opt_dynamics_lab[0])
#plot_evolution(opt_dynamics_lab)
In [ ]:
plot_cardinal([opt_dynamics_lab[0].states[0]])
plot_cardinal([opt_dynamics_lab[0].states[-1]])
Frequency analysis¶
In [ ]:
qubit_pulses = pulses_lab
time_list = tlist2
In [ ]:
qubit_pulses = opt_result.optimized_controls
time_list = tlist
In [68]:
def plot_spectrum(pulse, tlist, mark_freq=None, pos=1, xlim=None, mark_color=['k','k','k'], title="Pulse spectrum"):
samples = len(tlist)
sample_interval = tlist[-1]/samples
time = np.linspace(0, samples*sample_interval, samples)
signal_qubit = pulse
signal_spectrum = np.fft.fftshift(np.fft.fft(signal_qubit))
freqs = np.fft.fftshift(np.fft.fftfreq(samples, d=sample_interval))
plt.figure(figsize=(10,5))
plt.plot(freqs, np.abs(signal_spectrum)) # in GHz
if mark_freq is not None:
if not isinstance(mark_freq, list):
mark_freq = [mark_freq]
mf = np.array(mark_freq)/(2*π)
if pos==1:
plt.xlim(0, 2*mf[0])
elif pos==-1:
plt.xlim(-2*mf[0], 0)
elif xlim is not None:
plt.xlim(xlim[0]/(2*π), xlim[1]/(2*π))
[plt.axvline(x=m_f, ymin=0, ymax=1, color=col, linestyle='--', linewidth=1) for (m_f, col) in zip(mf,mark_color)]
plt.title(title)
plt.xlabel('f (GHz)');
plt.show();
plot_spectrum(qubit_pulses[0], time_list, mark_freq=[ω_q, ω_q + K_q, ω_q - K_q], pos=1)
plot_spectrum(qubit_pulses[1], time_list, mark_freq=[ω_q, ω_q + K_q, ω_q - K_q], pos=-1)